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Abstract — In recent work it is shown that importance 
sampling can be avoided in the particle filter through an 
innovation structure inspired by traditional nonlinear filtering 
combined with Mean-Field Game formalisms [9], [19]. The 
resulting feedback particle filter (FPF) offers significant variance 
improvements; in particular, the algorithm can be applied to 
systems that are not stable. The filter comes with an up-front 
computational cost to obtain the filter gain. This paper describes 
new representations and algorithms to compute the gain in the 
general multivariable setting. The main contributions are, 

(i) Theory surrounding the FPF is improved: Consistency 
is established in the multivariate setting, as well as well- 
posedness of the associated PDE to obtain the filter gain. 

(ii) The gain can be expressed as the gradient of a function, 
which is precisely the solution to Poisson's equation for 
a related MCMC diffusion (the Smoluchowski equation). 
This provides a bridge to MCMC as well as to approx- 
imate optimal filtering approaches such as TD-learning, 
which can in turn be used to approximate the gain. 

(iii) Motivated by a weak formulation of Poisson's equa- 
tion, a Galerkin finite-element algorithm is proposed for 
approximation of the gain. Its performance is illustrated 
in numerical experiments. 

I. Introduction 

In a recent work, we introduced a new feedback control- 
based formulation of the particle filter for the nonlinear 
filtering problem [18], [17]. The resulting filter is referred 
to as the feedback particle filter. In [18], [17], the filter 
was described for the scalar case, where the signal and the 
observation processes are both real-valued. The aim of this 
paper is to generalize the scalar results of our earlier papers 
to the multivariable filtering problem: 

AX, =a(X t )dt+dB t , (la) 
dZ, =h(X,)dt+dW t , (lb) 

where X, £ R d is the state at time f, Z, £ W" is the obser- 
vation process, a(-), h(-) are C 1 functions, and {B t }, {W,} 
are mutually independent Wiener processes of appropriate 
dimension. The covariance matrix of the observation noise 
{W t } is assumed to be positive definite. The function h is a 
column vector whose j-th coordinate is denoted as hi (i.e., 
h = (h\ ,!i2, . . . ,h m ) T ). For notational ease, the process noise 
{B t } is assumed to be a standard Wiener process. By scaling, 
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we may assume without loss of generality that the covariance 
matrices associated with {B,}, {W t } are identity matrices. 

The objective of the filtering problem is to estimate the 
posterior distribution of X, given the history % : = (j(Z s : s < 
t). The posterior is denoted by p*, so that for any measurable 
set A C R d , 

( p*{x,t)dx = ?{X t £A\%}. 
JxeA 

The filter is infinite-dimensional since it defines the evolu- 
tion, in the space of probability measures, of {p*(-,t) : t > 
0}. If a( • ), h( ■ ) are linear functions, the solution is given by 
the finite-dimensional Kalman filter. The article [3] surveys 
numerical methods to approximate the nonlinear filter. One 
approach described in this survey is particle filtering. 

The particle filter is a simulation-based algorithm to 
approximate the filtering task [6]. The key step is the 
construction of N stochastic processes {X t l : 1 < i < N}: The 
value X\ £ M. d is the state for the i th particle at time t. For each 
time t, the empirical distribution formed by, the "particle 
population" is used to approximate the posterior distribution. 
Recall that this is defined for any measurable set A C M. d by, 

P W (A,t) = ^l{Xl£A}. 

A common approach in particle filtering is called sequential 
importance sampling, where particles are generated accord- 
ing to their importance weight at every time step [1], [6]. 

In our earlier papers [18], [17], an alternative feedback 
control-based approach to the construction of a particle filter 
was introduced; see also [4], [14], [5], [12], [16] for related 
approaches. The resulting particle filter, referred to as the 
feedback particle filter, was described for the scalar filtering 
problem (where d = m = 1). The main result of this paper is 
to describe the feedback particle filter for the multivariable 
filtering problem ([Ta}-(|Tb"|>: 

The particle filter is a controlled system. The dynamics of 
the / th particle have the following gain feedback form, 

dX! =a(Xi)dt+dBi 

' W . '. (2) 

+ K{x; 7 t)dr,+Q(x;,t)dt, 

where {B' t } are mutually independent standard Wiener pro- 
cesses, /' is similar to the innovation process that appears in 
the nonlinear filter, 

dI l t :=dZ t - l -{h{xi)+h)dt, (3) 

where /; := E[h(XI)\3f t ]- In a numerical implementation, we 
approximate lf =1 h(Xj) =: U N l 



The gain function K is obtained as a solution to an 
Euler-Lagrange boundary value problem (E-L BVP): For 
j = l,2,...,m, the function <j)j is a solution to the second- 
order differential equation, 



V- (p(x,t)Vfy(x,t)) = -(hj(x)-hj)p(x,t) 



<j>j(x,t)p(x,t)dx = 0, 



(4) 



where p denotes the conditional distribution of X\ given 2£ t . 
In terms of these solutions, the gain function is given by 



(5) 



Note that the gain function needs to be obtained for each 
value of time t. 

Finally, ft = (Qi,ft2, —>£ld) is the Wong-Zakai correc- 
tion term: 



i d m 

ft/M^LLM^O^M- 



(6) 



The controlled system |2])-((6]l is called the multivariable 
feedback particle filter. 

The contributions of this paper are as follows: 

• Consistency. The feedback particle filter (|2]i is consistent 
with the nonlinear filter, given consistent initializations 
p(- ,0) = p*(-,0). Consequently, if the initial conditions 
{Xq}^ =1 are drawn from the initial distribution p*(-,0) of 
Xo, then, as N —> °°, the empirical distribution of the particle 
system approximates the posterior distribution p* ( • , t) for 
each t. 

• Well-posedness. A weak formulation of Q is introduced, 
and used to prove an existence-uniqueness result for (j)j in 
a suitable function space. Certain apriori bounds are derived 
for the gain function to show that the resulting control input 
in (2) is admissible (That is, the filter (2]) is well-posed in 
the Ito sense). 

• Numerical algorithms. Based on the weak formulation, a 
Galerkin finite-element algorithm is proposed for approxima- 
tion of the gain function K(x,t), The algorithm is completely 
adapted to data (That is, it does not require an explicit 
approximation of p(x,t) or computation of derivatives). 
Certain closed-form expressions for gain function are derived 
in certain special cases. The conclusions are illustrated with 
numerical examples. 

• Characterization of the feedback gain. The Smolu- 
chowski equation models a li-dimensional gradient flow with 
"noise": 

d* ( = -VSf(*,)df + V2d&, (7) 

where % is a standard Wiener process. It is regarded as 
the original MCMC algorithm: Under general conditions it 
is ergodic, with (unnormalized) stationary distribution '. 
The BVP (4) can be expressed as an instance of Poisson's 
equation for this diffusion, 



where @ is the differential generator for the Smoluchowski 
equation, with potential Sf = —log p. Subject to growth 
conditions on h and p, this implies the mean-integral rep- 
resentation for the vector-valued function, 



(j>j( X ) = J E[Ay(*,)-Ay|*o=Jt]df. 



(9) 



This representation also suggests an alternate proof of well- 
posedness and construction of numerical algorithms; cf., [8]. 
This will be the subject of future work. 

The outline of the remainder of this paper is as follows. 
The nonlinear filter is introduced and shown to be consistent 
in Sec|ll] The weak formulation of the BVP appears in Sec III 
where well-posedness results are also derived. Algorithms 
are discussed in Sec |IV] and a numerical example in Sec [V] 

II. Multivariable Feedback Particle Filter 



Consider the continuous time filtering problem (la lb I in- 
troduced in Sec U 

We denote as p*(x,t) the conditional distribution of X, 
given 3? t = a(Z s : s < t). The evolution of p*(x,t) is de- 
scribed by the Kushner-Stratonovich (K-S) equation: 



dp* = JzfV df + (h - h) T ( dZ, - h df )p* 



(10) 



@<j>j = -(hj-hj), l<j<m, 



(8) 



where h = J h(x)p*(x,t)6x and Jf^p* = —V • (pa) + \Ap, 
where A denotes the Laplacian in M d . 

A. Belief state dynamics & control architecture 
The model for the particle filter is given by, 

dX{ = a(Xi)dt+6Bl + u(Xl,t)dt + K(Xi,t)6Z t , (11) 

where X' t £ R d is the state for the 2 th particle at time f, 
and {B' t } are mutually independent standard Wiener pro- 
cesses. We assume the initial conditions {Xq}^ =1 are i.i.d., 
independent of {BJ}, and drawn from the initial distribution 
p*(x, 0) of Xo. Both {B l t } and {Xq} are also assumed to be 
independent of X t ,Z t . Note that the gain function K(x,t) is 
a d x m matrix and u(x,t) g R rf . 

We impose admissibility requirements on the control input 
Uj in {□}: 

Definition 1 (Admissible Input): The control input U\ is 
admissible if the random variables u(x,t) and K(x,t) 
are % = o(Z s : s < t) measurable for each t. And 
for each t, E[\u\] := E[£, h(X/,f)|] < °o and E[|K| 2 ] := 
E[L, ; -|Mx/,r)| 2 ]<~. ' ■ 

Recall that there are two types of conditional distributions 
of interest in our analysis: 

1) p(x,t): Defines the conditional dist. of X' t given 3? t . 

2) p*(x,t): Defines the conditional dist. of X t given 3%. 

The functions {u(x,t),K(x,t)} are said to be optimal if p = 
p*. That is, given p*(-, 0) = p(-,Q), our goal is to choose 
{m,K} in the feedback particle filter so that the evolution 
equations of these conditional distributions coincide (see ( [10] ) 
and (12)). 



X t = a{X t ) 
Vt = h(X t ) 



KI t 



Yt 



Y, 



■o- 



Xt = a(Xt) + Bi + K(Xi)li 
Y* = l(h(X}) + h) 



-o- 



(a): Kalman Filter (b): Feedback Particle Filter 

Fig. 1. Innovation error-based feedback structure for (a) Kalman filter and (b) nonlinear feedback particle filter (see Remark {TJ). 



The evolution equation for the belief state is described in 
the next result. The proof is identical to the proof in the 
scalar case (see Proposition 2 in [17]). It is omitted here. 

Proposition 1: Consider the process X' t that evolves ac- 
cording to the particle filter model (jTTJi. The conditional 
distribution of X} given the filtration p(x,t), satisfies the 
forward equation 

Ap = ^ pAt - V ■ (pK)AZ t 
- V • (pu) At 



d 

E 

l,k=\ 



dxid. 



iax k 



(MKK 7 



At. 



(12) 



B. General Form of the Feedback Particle Filter 

The general form of the feedback particle filter is obtained 
by choosing {m,K} as the solution to a certain E-L BVP 
based on p. The function K is a solution to 

V-(pK) = -(h-h) T p, (13) 

and the function u is obtained as 

u(x,t) = -^K(x,t) (h(x)+h) +Q(x,t). (14) 

The reader is referred to our earlier paper [17] for additional 
justification regarding these choices. 

Remark 1: Substituting (fl3|l-(fl"4"|i into ( fTT) gives the feed- 
back particle filter model <[2j-<[3j in Sec [I] 

In the Stratonovich form, the filter admits a simpler 
representation, 



dX! =a(x;)At + AB' t + K(X',t)o ( dZ, - -(A (3?) +h)At 



Given that the Stratonovich form provides a mathematical 
interpretation of the (formal) ODE model [15, Section 3.3 
of the SDE text by 0ksendal], we also obtain the (formal) 
ODE model of the filter. Denoting Y t 

~d7 



^ and white noise 
process B\ = ^j-, the ODE model of the filter is given by, 



AX[ 
At 



a(X t ') +B\ + K(X',t) [Y t - -(h(x;)+h) 



1 



The feedback particle filter thus provides a generalization of 
the Kalman filter to nonlinear systems, where the innovation 
error-based feedback structure of the control is preserved (see 
Fig. [T]). For the linear case, it is shown in Sec III-D that the 
gain function is the Kalman gain. For the nonlinear case, 
the Kalman gain is replaced by a nonlinear function of the 
state. ■ 



If one further assumes that the control input U' t is admissi- 
ble, a short calculation shows that the feedback particle filter 
is consistent with the choice of {m,K} given by (fT~3l> - (fl~4]> . 
This calculation appears in Appendix [A] 

C. Consistency with the Nonlinear Filter 

To establish admissibility of the input U\ requires addi- 
tional assumptions on the density p and function h: 

(i) Assumption Al The probability density p(x,t) is of 
the form p(x,t) = e~^( x >'\ where &(x,t) is a twice 
continuously differentiable function with 

|VSf| 2 -2ASf as|x|^oc. (15) 

(ii) Assumption A2 The function h satisfies, 



\h(x)\ 2 p(x,t)Ax < oo. 



where |/j(-x)| 2 := \hj(x)\ 2 . For admissibility of u, our 
arguments require additional assumptions: 
(iii) Assumption A3 The second derivatives of @(x,t) 
with respect to x are uniformly bounded at each t, i.e., 
for all x € M. d , t > 0. 



xjdx k 



] g^(x,t)\<C 2 (t) 



(iv) Assumption A4 The first (weak) derivatives of h 
satisfy 

\Vh{x)\ 2 p{x,t)Ax<™, 

where |V^)| 2 := ^ |g(x)| 2 . ■ 

Under these assumptions, it is shown in Theorem[2]that the 
gradient-form Q of the E-L BVP ( fT3| is uniquely obtained 
to give and thence K. 

The admissibility of the resulting control input is estab- 
lished in Corollary JT1 Theorem [2] and Corollary [T] are stated 
and proved in Sec [Hi] 

The following theorem then shows that the two evolution 
equations ( fT0) i and (jT2]» are identical. The proof appears in 
Appendix [A] 

Theorem 1: Consider the two evolution equations for p 
and p*, defined according to the solution of the forward equa- 
tion ( fT2l > and the K-S equation ( fT0| >, respectively. Suppose 
that the gain function K(x,f) is obtained according to (|4|-(|5]l. 
Then, provided p( ■ ,0) = p*( ■ ,0), we have for all t > 0, 

p(-,t)=p*(-,t). 



III. Existence, Uniqueness and Admissibility 

The aim of this section is to introduce a particular 
gradient-form solution of the BVP (13) . The gradient-form 
solution is obtained in terms of m real-valued functions 
{fa(-,t),<h(-,t),...,bn(-,t)}. For ;= 1,2,. ..,m, the func- 
tion is a solution to, 



V-(p(x,t)V<j) j (x,t)) = -(hj(x)-tij)p(x,t), 



(j)j(x,t)p(x,t)dx = 0. 



(16) 



The normalization / <j>j(x,t)p{x,t)dx = is for convenience: 
If (j>j is an solution to the differential equation (To} , we obtain 
the desired normalization on subtracting its mean. 

In terms of these solutions, the gain function is given by, 



Kij(x,t) 



dxi 



(x,t), 



x£ 



(17) 



It is straightforward to verify that K thus defined is a 
particular solution of the BVP (13) . 

A. Poisson's Equation Interpretation 

The differential equation ( [To} is solved for each t to give 
the m functions {(j>j{-,t) : 1 < j < m}. On dividing each 
side of this equation by p, elementary calculus leads to the 
equivalent equation ([SJ, with generator @ defined for C 2 
functions / via, 

0/=-V#-V/+A/ 

and with Sf(-) = — logp(- ,t). This is the differential gener- 
ator for the Smoluchowski equation Q. It is shown in [10] 
that this diffusion is exponentially ergodic under mild condi- 
tions on Sf. Consequently, E[hj(<P t ) — hj | <t>o =x] converges 
to zero exponentially fast, subject to growth conditions on 
hj, and from this we can conclude that |9]l is well defined, 
and provides a solution to Poisson's equation ([8]i [8]. 

Poisson's equation can be regarded as the value function 
that arises in average-cost optimal control, and this is the 
object of interest in the approximation techniques used 
in TD-learning for average-cost optimal control [13]. The 
integral representation |9} suggests approximation techniques 
based on approximate models for the diffusion <f>. 

The remainder of this section is devoted to showing 
existence and uniqueness of the solution of ( [To} , and ad- 
missibility of the resulting control input, obtained using gain 
function defined by (17) . 

B. Weak Formulation 

Further analysis of this problem requires introduction of 
Hilbert spaces: L 2 (R d ;p) is used to denote the Hilbert space 
of functions on M. d that are square-integrable with respect 
to density p(-,t) (for a fixed time f); H k (R d ;p) is used to 
denote the Hilbert space of functions whose first /c-derivatives 
(defined in the weak sense) are in L 2 (Mr;p). Denote 



//(HIT;/?) := <^ ^eH l (W;p) / <j)(x)p(x,t)dx = 



A function 0; G Hq (R d ; p) is said to be a weak solution 
of the BVP (0 if 

/ V(j>j(x,t) ■V\i/(x)p(x,t)dx = / (hj(x) — fij)y(x)p(x,t)dx, 

(18) 

for all y eH l (R d ;p). 

Denoting E[] := / -p(x,t)dx, the weak form of the 
BVP (To} can also be expressed as 

E[V<l>j-V\ir] = E[(hj-hj)\ir], V y/ e H 1 (M d ;p). (19) 

This representation is useful for the numerical algorithm 
described in Sec HVl 

C. Main Results 

The existence-uniqueness result for the BVP (To} is de- 
scribed next — Its proof is given in Appendix [B] 

Theorem 2: Under Assumptions A1-A2, the BVP (To} 
possesses a unique weak solution (j)j G Hq (Mr; p), satisfying 



|V^-(x)| 2 ;>(jc,f)dx < T / \hj(x)-hj\ 2 p(x,t)6x. (20) 



If in addition Assumptions A3-A4 hold, then <j>j G H 2 (R d ; p) 
with 

J \(D 2 <l)j)(V(j)j)\p(x 1 t)dx<C(X;p) J \Vhj\ 2 p(x,t)dx, 

(21) 

where X is (spectral gap) constant (see Appendix [B]i and 

C(A;p) = ^(^*+l)'' 



1/2 



a 3 / 2 v x 

The apriori bounds (20}-(2T} are used to show that the 
control input for the feedback particle filter is admissible. 
The proof is omitted on account of space. 

Corollary 1: Suppose <j>j is the weak solution of BVP (To} 
as described in Theorem [2] The gain function K is obtained 
using (TT} and u is given by (14) . Then 

i m 

E[|K| 2 ]< r £ 

.7=1 
1 

A 4 



E[|h|1< 



\hj(x)\ 2 p(x,t)dx, 



\Vhj\ 2 )p(x,t)dx, 



where C(X;p) is given in Theorem |5] That is, the resulting 
control input in (T7} is admissible. ■ 



D. Linear Gaussian case 
Consider the linear system, 

dX, =AX,dt+dB, 
dZ,=HX,dt+dW t 



(22a) 
(22b) 



where A is an d x d matrix, and H is an m x d matrix. The 
initial distribution p* (x, 0) is Gaussian with mean vector jio 
and co variance matrix £o- 

The following proposition shows that the Kalman gain is 
a gradient-form solution of the multivariable BVP (T3j ): 

Proposition 2: Consider the d-dimensional linear sys- 
tem (22a[ )- (22b} . Suppose p(x,t) is assumed to be Gaussian: 
p(x,t) = A i expi-Ux- pi t ) T Y.y\x- ji t )), where x = 

(27r)2|I,| 3 



(x u x 2 , ■■■,x d ) T , fi, = (jUit,M2n ■ ■ -,IJ-dt) T is the mean, L, is the 
covariance matrix, and |E,| > denotes the determinant. A 
solution of the BVP ( |T6] > is given by, 

d 

fyfoO = Y J %H T ] kj {x k -iik,)- (23) 

Using ( |17) , K(x,t) = L t H T (the Kalman gain) is the gradient 
form solution of (jT3]l. ■ 

The formula ( |23| l is verified by direct substitution in the 
BVP (jT6]» where the distribution p is multivariable Gaussian. 

The gain function yields the following form for the particle 
filter in this linear Gaussian model: 

dX; = A Xj dt + dB\ + T.,H T ( dZ, - H *' ^ At \ . (24) 

Now we show that p = p* in this case. That is, the con- 
ditional distributions of X and X' coincide, and are defined 
by the well-known dynamic equations that characterize the 
mean and the variance of the continuous-time Kalman filter. 

Theorem 3: Consider the linear Gaussian filtering prob- 
lem defined by the state-observation equations ( |22a[ )-( |22b| . 
In this case the posterior distributions of X t and X\ are 
Gaussian, whose conditional mean and covariance are given 
by the respective SDE and the ODE, 

dju, = AjU, df + T. t H T ( dZ, - H/i, dt) 

— E, = AT., + T.,A T + 1 - T.,H T HI. t 
dt 

m 

The result is verified by substituting p(x,t) — 
(27t)~z|£,|~iexp [-\(x- pi,) T Y.^ l (x- pi,)] in the forward 
equation ( [12) . The details are omitted on account of space, 
and because the result is a special case of Theorem [T] 

In practice {jj, t ,T. t } are approximated as sample means and 
sample covariances from the ensemble {X}}f =1 : 

1 N 

u ~ u {N] — -VX' 
The resulting equation |24| for the ; th particle is given by 

dx; = a x; dt + dB\ + Y.f ] H T dz t - h ' 2 Mf dt . 

As — > oo, the empirical distribution of the particle system 
approximates the posterior distribution p*(x,t) (by Theo- 
rem [3). 

IV. Finite-Element Algorithm 

In this section, a Galerkin finite-element algorithm is 
described to construct an approximate solution of ( |18) . Since 
there are m uncoupled BVPs, without loss of generality, we 
assume scalar-valued observation in this section, with m=l, 
so that K = V(j>. The time t is fixed. The explicit dependence 
on time is suppressed for notational ease (That is, p(x,t) is 
denoted as p(x), $(x,t) as 0(x) etc.). 



A. Galerkin Approximation 

Using ( fT9| ), the gain function K = V<p is a weak solution 

if 

E[K-Vy] = E[(h-h)y], V y e H l (R d ; p) . (25) 
The gain function is approximated as, 

L 

K= Y,KiXi(x), 

i=i 

where {Xi( x )}f=i are basis functions. For each / = 1,...,L, 
Xi(x) is a gradient function; That is, %i{x) = VC,i(x) for some 
function Q(x) G H^{R d ;p). 

The test functions are denoted as {Wk(x)} k=i and S := 
span{y/i (x), y/ 2 (x), . . . , y/ L (x)} c H 1 (R d ;p). 

The finite-dimensional approximation of the BVP ( |25| ) is 
to choose constants {K/}f =l such that 

L 

£jc / Eto-VvA] = E[(A-A)vA] ! V v/ e 5. (26) 

i=i 

Denoting \A] kl = ■ V%], b k = E[(h - tyy k ], K = 
{k\ , JC2, . . . , Kl) t , the finite-dimensional approximation (|26j> 
is expressed as a linear matrix equation: 

Ak = b. 

The matrix A and vector b are easily approximated by using 
only the particles: 

[A] k , = E\Xi ■ V % ] w i £ • VW^), (27) 

ft 4 = E[(h-h) Wk ] » ^£(A(X/)-A)%(X/), (28) 
iV i=i 

where recall ft » ^ J^j ft(X/'). 

5. Example 1: Constant Gain Approximation 

Suppose Xi — e h the canonical coordinate vector with 
value 1 for the / th coordinate and zero otherwise. The test 
functions are the coordinate functions y/ k (x) = x k for k = 
1,2,. ..,d. Denoting y/(x) = (y/i , 1/2, ■ ■ • , Wd) T = x, 

K = E[K] = E[(h-h)Y] = J (h(x)-h)y(x)p(x)dx 

*ltm)-h)Xl- (29) 
ly i=i 

This formula yields the constant-gain approximation of the 
gain function. 

C. Example 2: Single-state Case 

Consider a scalar example, where the density is a sum of 
Gaussian, 

3 

p(x) w £AVM, 

where qj(x) = q{x;^^) = 7=jexp(-M£), V > 0, 
£A' = 1. The parameter values for A^/i^E' are tabulated 
in Table I. 



(a) Table 1 
Parameter Values 



3 


X j 




E j 


1 


0.2 


-1 


0.25 


2 


0.5 





0.25 


3 


0.3 


1 


0.25 



0.5 



- (b) 








_,,-'''''lpl (x) 


30 






20 







0.5 



Fig. 2. (a) Parameter values and (b) (V1.V2) m the Example. 

In the scalar case, a direct numerical solution (DNS) of 
the gain function is obtained by numerically approximating 
the integral 

1 



K(x) = -— / (h(y)-h)p(y)dy. 

P(X) J -00 

The DNS solution is used to provide comparisons with the 
approximate Galerkin solutions. 

The Galerkin approximation of the gain function is con- 
structed on an interval domain Del. The domain is a union 
of finitely many non-intersecting intervals D/ — [a/_i,a/), 
where oq < a\ < ... <ul- 

Define for / = 1,2, . . . ,L and k = 1,2, ... ,L: 

Basis functions: Xi( x ) — ^D t {x), 
Test functions: Wk( x ) = \x — ak\- 

Figure |2] depicts the test functions {y/i (jc), y/it*)} for D = 
[0, 1] and ao = 0, a\ = j and a 2 = 0. The basis functions are 
the indicator functions on [0, 5) and [j,l). 

Figure [3] depicts a comparison of the DNS solution and the 
Galerkin solution for h(x) ~x 2 , D = [—2,2] and L = 1,5, 15. 
For a given L, the basis and test functions are constructed for 
a uniform partition of the domain (That is, a/ = —2 + ^4). 
The Galerkin solution is obtained using N = 1000 particles 
that are sampled from the distribution p. The particles are 
used to compute matrix A and vector b, using formulae p7| ) 
and ( |28) , respectively. Since the analytical form of p is 
known, these matrices can also be assembled by using the 
integrals: 



\Ahi = J Xi{x)-^Yk{x)p(x)dx, 
bk = / (h(x) — h)yk(x)p(x)dx. 



(30) 
(31) 



The figure also depicts the Galerkin solution based on the 
integral evaluation of the matrix A and vector b. 

For L — 15, the matrix A was found to be singular for the 
particle-based implementation. This is because there are no 
particles in D15. In this case, the Galerkin solution is obtained 
using only the integral formulae (|30l)-(|3T|). These formulae 
are exact while the particle-based formulae ( p7] i and ( f28| are 
approximations. In the other two cases (L = 1 and L = 5), 
the particle-based solution provides a good approximation. 

V. Numerics 

Consider a target tracking problem with two bearing-only 
sensors [2]. A single target moves in a two-dimensional 
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Fig. 4. Simulation results: Comparison of the true target trajectory with 
the estimate obtained using FPF 



(2d) plane according to the standard white-noise acceleration 
model: 

AX, =AX t At + TAB u 

where X := (X u VuX 2 ,V 2 ) T G M 4 , (X U X 2 ) denotes the posi- 
tion and (Vi,V2) denotes the velocity. The matrices, 



10 



1 







1 

1 



and B t is a standard 2d Wiener process. 
The observation model is given by, 

dZ, =h{X,)At + a w AW t , 

where W, is a standard 2d Wiener process, h — {hi^h^f and 



hj(x\, vi,x 2 , V2) = arctan 



x 2 X2 



sen j) * 



X\ —x 



sen j) 



j=l,2, 



where (x 



( sen ;) „( sen /h 



j A 2 



denote the position of sensor j. 



Figure [4] depicts a sample path obtained for a typical 
numerical experiment. The sensor and target locations are 
depicted together with an estimate (conditional mean) that 
is approximated using a feedback particle filter. The back- 
ground depicts the ensemble of observations that were made 
over the simulation run. Each point in the ensemble is 
obtained by using the process of triangulation based on 
two (noisy) angle measurements. The simulation parameters 
are: The initial position of the target is depicted, the initial 
velocity was chosen as (0.2,-5) and Ob = 0.1; The two 
sensor positions are depicted and c% = 0.017; The particle 
filter comprised of N — 200 particles whose initial position 
was chosen from a Gaussian distribution whose mean is 
depicted. The gain function was obtained using the constant 
gain approximation in ( |29) , The simulation results show that 
the filter can adequately track the target. 

Appendix 

A. Proof of Theorem [7] 

It is only necessary to show that with the choice of {u, K} 
given by (p~3])-(|T4]>^ we nave d /?(x,f) = Ap*(x,t), for all x and 




f, in the sense that they are defined by identical stochastic 
differential equations. Recall dp* is defined according to 
the K-S equation ( fTO] ), and dp according to the forward 
equation (12) . 

Recall that the gain function K is a solution of the 
following BVP: 

V-(pK) = -p(h-h) T . (32) 
On multiplying both sides of ( fT4| i by — p, we obtain 

—up — ijj>{h — /i)K — Q.p + pKh 

= -^K[V-(/?K)] r -Qp+pKh 

where ( |32] > is used to obtain the second equality. Denoting 
E := jK[V- (pK)] T , a direct calculation shows that 

E l +n l p= J E^-{p[KK T ] lk ). 

k=\ " Xk 



(33) 



Substituting this in ( |33| >, on taking the divergence of both 
sides, we obtain 

1 



.. 5 -(p[KK r ] tt )=V-(pK)ft. (34) 

Using < [32] > and p4| in the forward equation ( |12) , 

dp = jSf f /> + (/i - ( dZ, - A d/)p. 
This is precisely the SDE ( fTO) , as desired. ■ 

B. Proof of Theorem [2] 

We omit the subscript "/' in this proof, writing just (j) and 
h. We also suppress explicit dependence on time t, writing 
p{x) instead of p(x,t) and <j)(x) instead of $(x,t). 

Under Assumption Al, p is known to admit a spectral gap 
(or Poincare inequality) with constant A > ci [11]: That is, 
for all functions (j) G H^W 1 ; p), 

j \<j>(x)\ 2 p(x)dx< i j \V<^(x)\ 2 p{x)dx. (35) 
Consider now the inner product 

<</>,y/>:= V(j)(x) ■V\j/(x)p(x)dx. 



On account of ( |35| ), the norm defined by using the inner prod- 
uct < •, • > is equivalent to the standard norm in H l (R d ;p). 



(i) Consider the BVP in its weak form ( fT8j ). The integral 
on the righthand-side is a bounded linear functional on 
y/ € Hq, since 

— h)y/(x)p(x)dx\ 2 
< f \h{x)-h\ 2 p{x)dx j \\\r(x)\ 2 p{x)dx 



< (const.) / \Vy/(x)\ p(x)dx, 



where ( f33p is used to obtain the second inequality. 
It follows from the Reisz representation theorem that 
there exists a unique <p € Hq such that 

<</>,y/>= j(h(x)—h)y(x)p(x)dx, 

for all v/ € //J p). Thus is a weak solution of the 
BVP, satisfying ( fTg) . 

(ii) Suppose is a weak solution. Using y/ = <p in ( fT8| , 



|V0| 2 p(jt)dx = J (h(x)-h)<j)(x)p(x)dx 
<[ [ \h(x)-h\ 2 p(x)dx Y ( [ \<$>{x)\ 2 p{x)dx 



[ j \hl.v)-fi\ 2 p(x)dx) 2 (j / |V^(jc)| 2 p(.vidv 



by f35| >. The estimate (|20J follows, 
(iii) For the final estimate ( f2T| , we need: 
Lemma 1: Under Assumptions A1-A4, the weak solu- 
tion of the BVP ([18]) belongs to H 2 (R d ;p), with 



/ 



\D 2 (j>\ 2 pdx< j V^-Gpdx (36) 

where the vector function G G L 2 (M'';/?) is defined by 
G = D 2 (logp)V^+V/! 

and where \D 2 <j>\ 2 = Zj, ki^f- 

Proof: First note that each entry of the Hessian 
matrix D 2 (logp) is bounded, by Assumption A3, and 
that V/i e L 2 (R d ;p) by Assumption A4. Hence G G 
L 2 (R;/j). 

Next, elliptic regularity theory [7, Section 6.3 of the 
PDE text by Evans] applied to the weak solution G 



says that e Hf 



differential equation holds pointwise: 

-V-(/7V0) = (h-h)p. 
We differentiate with respect to x k to obtain: 

dx k 



. Hence the partial 
(37) 



■ V.(pV^)-V( 



5l0g ^-W)-^v. W ) 



dx k 



^—p+{h-h)— p. 

dx k dx k 



The final terms on the left and right sides cancel, by 
equation < [37] >. Thus the preceding formula becomes 



-V.(pV 



50 

5xt 



GkP, 



(38) 



where denotes the kth component of G(x). 
Let (x) > be a smooth, compactly supported "bump" 
function, meaning j8(x) is radially decreasing with 
j8(0) = 1. Let s > and multiply ((38) by j3(m) 2 ^. 
Integrate by parts on the left side (noting the boundary 
terms vanish because j3 has compact support) to obtain 



/ 



dxk 



G k pdx. 



The right hand side RHS — > J ^-G k pdx by dominated 
convergence as s->0, since j3 (0) = 1 . The left side, 



LHS = / J3(«) 2 |v|^| 2 pdx 



dtp 



dx k 



■2s I ^P(sx)(VP)(sx) 
dx k 



dx k 



Clearly the second term is bounded by 

2*||Vj3|U E 



^\P(sx)\^\ P 6x 
dx k dx k 



and so 

(l-s\\Vp\\ L ~ (v 



dx k 



-P(sx) 2 \V 



dx k 



pdx 



50, 
dx k 



P{sx) 2 \V^\ 2 pdx 
dx k 

2 dx<LHS. 



Letting s -> in LHS and RHS, and recalling that /5(jc) 
is radially decreasing, we conclude from the monotone 
convergence theorem that 

dtp 



|V^| 2 ,dx< 
dx k 



1X k 



G k pdx. 



Summing over k completes the proof of the lemma. 
Next we prove pT) . We will use several times that 



/ 



|V0| 2 pdx< -j \h-h\ 2 pdx<-^ 



\Vh\ 2 pdx, 



(39) 

by ( |20| ) followed by ( (35) applied to the function h—h€ 



We have 



{D 2 (j))(V(j))\pdx< J |D 2 0||V0|pdx 

<(l 2 J\Vh\ 2 pdx)y\j\G\ 2 pdx)^ 

by Lemma [T] Cauchy-Schwarz, and ( |39) . The defini- 
tion of G, the L 2 -triangle inequality and ((39) show that 



(/|G| 2 pd*) 



1/2 

^l^aog^lkC/iV^P/^dv) 1 2 H-(/|V/,i 2 /Hlvi' 2 
^ 2 (logp)||i 



< 



1 ( / \Vhfpdx) 



1/2 



The estimate ( |2T) now follows. 
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